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I.  INTRODUCTION 

Since  the  very  beginning  of  the  application  of  numerical  methods  to  elec¬ 
tromagnetics,  there  has  been  a  keen  interest  in  developing  computer  codes  for 
treating  radiation  and  scattering  problems  involving  arbitrarily-shaded  con¬ 
ducting  bodies.  Of  the  various  possible  approaches  available  to  developers 
of  such  codes,  the  most  commonly  used  have  been  wire-grid  and  surface  patch 
modeling  in  conjunction  with  integral  equation  formulations. 

The  wire-grid  modeling  approach  has  been  remarkably  successful  in  many 
problems,  particularly  in  those  requiring  the  prediction  of  far-field  quanti¬ 
ties  such  as  radiation  patterns  and  radar  cross-sections  [1].  The  approach  is 
not  as  well-suited  for  calculating  near-field  and  surface  quantities,  such  as 
surface  current  and  input  impedance,  however.  Some  of  the  difficulties  en¬ 
counter  i  in  wire-grid  modeling  include  the  occasional  ■presence  of  ficticious 
loop  currents  in  the  solution,  difficulties  with  internal  resonances  [2],  and 
problems  of  relating  computed  wire  currents  to  equivalent  surface  currents. 

The  accuracy  of  wire-grid  modeling  has  also  been  questioned  on  theoretical 
grounds  [3].  These  difficulties  have  provided  strong  incentives  for  developing 
surface  patch  approaches  as  alternatives  to  wire-grid  techniques. 

Several  approaches  to  surface  patch  modeling  have  been  suggested.  Knepp 
and  Goldhirsh  [4]  partition  a  surface  into  non-planar  quadrilateral  patches 
and  employ  the  magnetic  field  integral  equation  (MFIE)  to  solve  the  electro¬ 
magnetic  problem.  Albertsen  et  al.  [5]  solve  for  the  current  and  compute  ra¬ 
diation  patterns  for  satellite  structures  with  attached  wire  antennas,  booms, 
and  solar  panels .  They  use  the  MFIE  with  planar  quadrilateral  surface  patches 
to  model  the  satellite,  and  use  the  electric  field  Integral  equation  (EFIE) 
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to  treat  the  wire  antennas  in  their  hybrid  formulation.  The  arbitrary  surface 
treatment  of  the  widely-used  Numerical  Electromagnetic  Code  (NEC)  developed 
at  the  Lawrence  Livermore  Laboratory  [6]  is  also  based  on  the  formulation  of 
Albertsen  et  al.  Wang  et  al.  [7]  extend  the  use  of  piecewise-sinusoidal  basis 
functions,  well-known  in  wire  analyses,  to  the  treatment  of  surfaces.  They 
use  an  EFIE  formulation  and  model  surfaces  by  planar  rectangular  patches. 

Sankar  and  Tong  [8]  employ  planar  triangular  patches  to  model  a  square  plate 
and  point  out  the  possibility  of  extending  their  approach  to  arbitrary  bodies. 
Their  formulation,  based  on  a  variational  formula  for  the  current  which  is 
made  stationary  with  respect  to  a  set  of  trial  functions,  is  equivalent  to  a 
Galerkin  solution  of  the  EFIE.  Wang  [9,10]  employs  planar  triangular  patches 
in  conjunction  with  the  MFIE  and  uses  basis  functions  which  contain  the  phase 
variation  of  the  incident  field  in  each  jatch.  Unfortunately,  this  procedure 
makes  the  resultant  moment  matrix  depend  on  the  incident  field.  Jeng  et  al. 
[11]  propose  using  the  MFIE  and  non-planar  triangles  to  model  arbitrary  sur¬ 
faces.  Singh  and  Adams  [12]  propose  using  planar  quadrilateral  patches  and 
sinusoidal  basis  functions  with  the  EFIE. 

In  arbitrary  surface  modeling,  the  EFIE  has  the  advantage  that  It  applies 
to  both  open  and  closed  bodies,  whereas  the  MFIE  applies  to  closed  bodies  only. 
On  the  other  hand,  for  arbitrarily-shaped  objects  the  EFIE  is  much  more  diffi¬ 
cult  to  deal  with,  as  attested  to  by  the  fact  that  of  the  EFIE  formulations 
discussed,  only  Wang  et  al.  have  actually  treated  non-planar  structures — and 
their  formulation  is  limited  to  structures  with  curvature  in  one  dimension 
only.  The  difficulties  with  the  EFIE  stem  primarily  from  the  presence  of 
derivatives  and  a  singular  kernel  in  the  integral  equation.  One  manifestation 
of  the  derivatives  is  that  if  basis  functions  representing  the  current  are  not 
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constructed  such  that  their  normal  components  are  continuous  across  surface 
edges,  then  line  or  point  charges  are  deposited  along  the  edges.  If  present, 
these  ficticious  charges  usually  lead  to  deleterious  effects  in  the  solution. 

The  approach  of  Wang  et  al.  [7]  is,  as  they  point  out,  free  of  these  diffi¬ 
culties,  but  their  use  of  rectangular  patches  restricts  their  consideration 
to  (finite)  cylindrical  surfaces.  More  appropriate  for  modeling  arbitrarily- 
shaped  surfaces  are  planar  triangular  patch  models  such  as  shown  in  Fig.  1. 

Some  of  the  advantages  of  triangular  patch  surface  modeling  have  been 
noted  by  Sankar  and  Tong  [8],  as  well  as  Wang  [9].  For  example,  triangular 
patches  have  the  ability  to  conform  to  any  geometrical  surface  and  boundary, 
they  permit  simple  descriptions  of  the  surface  and  parch  scheme  to  the  computer, 
and  they  may  be  used  with  greater  patch  densities  on  those  portions  of  the 
surface  where  more  resolution  is  desired.  (Although  planar  quadrilateral 
patch  modeling  shares  many  of  these  features,  the  vertices  of  planar  quadri¬ 
laterals  may  not  be  independently  specified  because  all  four  vertices  must 
lie  in  the  same  plane.) 

In  this  report  we  use  planar  triangular  patch  modeling  and  apply  the 
method  of  moments  [13]  to  develop  numerical  procedures  for  both  the  EFIE  and 
MFIE  formulations.  The  computer  code  based  on  the  EFIE  is  capable  of  handling 
either  open  or  closed  and  arbitrarily-curved  structures  of  finite  extent. 
Discounting  limitations  of  the  computer,  the  code  can,  in  fact,  treat  any 
object  whose  surface  is  orientable,  connected  (i.e.,  the  body  does  not  com¬ 
prise,  in  reality,  two  or  more  separate  objects),  and  free  of  intersecting 
surfaces.  Not  only  open  and  closed  surfaces,  but  also  multiply-connected 
objects  such  as  the  structure  with  a  "handle"  (c.f.  Appendix  A)  shown  in  Fig. 

1  are  admissable.  The  computer  code  based  on  the  MFIE  has  the  same  range  of 
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Fig.  1.  Arbitrary  surface  modeled  by  triangular  patches. 
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applicability,  except  that  it  is  restricted  to  closed  surfaces.  Both  the 
EFIE  and  MFIE  approaches  developed  are  simple  and  efficient  to  apply. 

We  remark  that  a  previous  paper  has  considered  the  electrostatic  problem 
of  determining  the  charge  distribution  on  arbitrarily-shaped  conducting  bodies 
modeled  by  triangular  patches  [14].  We  note  also  that  the  formulation  used  there 
is  related  to  the  static  limit  of  the  present  EFIE  formulation  and  that  both 
formulations  employ  piecewise  constant  charge  representations. 

In  the  following  section,  we  present  the  EFIE  formulation.  A  new  set  of^ 
basis  functions  defined  on  triangular  patches  is  described  and  used  to  repre¬ 
sent  the  current  in  the  moment  method.  These  new  basis  functions  are  free  of 
ficticious  line  or  point  charges  and  are  analogous  to  the  so-called  "rooftop" 
functions  defined  on  planar  rectangular  subdomains  [15]. 

The  MFIE  formulation,  which  also  makes  use  of  these  new  basis  functions, 
is  presented  in  Section  III.  Despite  its  lack  of  generality,  the  MFIE  formu¬ 
lation  is  important  because  its  moment  matrix  is  required  in  problems  of 
scattering  by  dielectric  objects  [16],  and  in  the  so-called  combined  field 
integral  equation  formulation  [17].  The  latter  is  a  technique  for  eliminating 
difficulties  in  both  the  EFIE  and  MFIE  formulations  for  scattering  problems 
at  frequencies  corresponding  to  the  cavity  resonances  of  the  interior  region 
for  closed  surfaces. 

In  Section  IV,  numerical  results  obtained  using  the  El'IE  formulation  are 
presented  for  triangular  patch  models  of  a  flat  square  plate,  a  bent  rectang¬ 
ular  plate,  a  circular  disk,  and  a  sphere.  Results  obtained  using  the  MFIE 
are  also  presented  for  the  sphere  problem. 
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II.  ELECTRIC  FIELD  FORMULATION 

In  this  section,  we  derive  an  integral  equation  for  the  surface  current 
induced  on  a  conducting  scatterer  from  the  boundary  conditions  on  the  electric 
field.  A  set  of  expansion  functions  and  a  testing  procedure  are  then  devel¬ 
oped  for  use  in  applying  the  method  of  moments,  and  the  moment  matrix  is  de¬ 
rived.  Finally,  the  evaluation  of  elements  of  the  moment  matrix  is  discussed. 

Electric  Field  Integral  Equation 

Let  S  denote  the  surface  of  an  open  or  closed,  perfectly  conducting 
scatterer.  An  electric  field  E1,  defined  in  the  absence  of  the  scatterer, 
is  incident  and  induces  surface  currents  J  on  S.  If  S  is  open,  we  regard  J 
at  each  point  as  che  vector  sum  of  the  currents  on  opposite  sides  of  the  sur- 

—  s 

face.  We  can  compute  the  scattered  electric  field  E  from  the  surface  current 
by 


ES=  -jioA  -  V41 

where  the  magnetic  vector  potential  is  defined  as 

u  f  _  e~jkR 

=  4¥  J  J  V-  ds’ 

S 

and  the  scalar  potential  is 

1  f  e"jkR 

*<r>  ‘  4^J  °  ~R“  dS'- 
s 


(1) 


(2) 


(3) 
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An  exp(jwt)  time  dependence  is  assumed  and  suppressed,  and  k  =  =  2n/X, 

where  A  is  the  wavelength.  The  permeability  and  permittivity  of  the  sur¬ 
rounding  medium  are  y  and  e,  respectively,  and  R  =  |r  -  r'|  is  the  distance 
between  an  arbitrarily-located  observation  point  r  and  a  source  point  r'  lo¬ 
cated  on  S.  Both  r  and  r'  are  defined  with  respect  to  a  global  coordinate 
origin  0.  The  surface  charge  density  0  is  related  to  the  surface  divergence 
of  J  through  the  equation  of  continuity, 

Vg-  J  =  — j 0X7.  (4) 

We  derive  the  integro-dif f erential  equation  for  J  by  applying  the  boun- 
~ ”  i  ”  s  — 

dary  condition  n  x  (E  +  E  )  =  0  on  S,  obtaining, 

-EL  -  W^-Wtan-  1  S-  (5> 

Eq.  (5),  with  (2)  -  (4),  constitutes  the  so-called  electric  field  integral 
equation  (EFIE).  One  notes  that  the  presence  of  derivatives  on  the  current 
in  (4)  and  on  the  scalar  potential  in  (5)  suggests  that  one  should  be  careful 
in  selecting  the  expansion  functions  and  testing  procedure  in  the  method  of 
moments.  In  the  next  section,  we  choose  expansion  functions  which  yield  a 
continuous  current  and  a  piecewise  constant  charge  representation. 

Deve lopment  o f  Bas is  Functions 

In  this  section,  we  discuss  a  set  of  basis  functions,  originally  pro¬ 
posed  by  Glisson  [18],  which  are  suitable  for  use  with  the  EFIE  and  triang¬ 
ular  patch  modeling.  We  assume  a  suitable  triangulation  approximating  S  and 
defined  by  a  set  of  faces,  edges,  vertices,  and  boundary  edges  (c.f.  Fig.  1) . 
Fig.  2  shows  two  triangles,  T^  and  T^,  associated  with  the  n  edge  of  a 
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triangulated  surface  modeling  a  scatterer.  Points  in  may  be  designated 

either  by  the  position  vector  r,  defined  with  respect  to  0,  or  by  the  position 

vector  p+,  defined  with  respect  to  the  free  vertex  of  T+.  Similar  remarks 
n  n 

apply  to  the  position  vector  except  that  it  is  directed  toward  the  free 
vertex  of  T^.  It  is  assumed  that  the  plus  or  minus  designation  of  the  tri¬ 
angles  has  been  chosen  such  that  the  positive  current  reference  direction 

(c.f.  Appendix  A)  associated  with  the  n*"*1  edge  is  from  T+  to  T  .  We  define 

n  n 

a  vector  basis  function  associated  with  the  n^  edge  as 


fn(?)  = 


i  + 

where  Z  is  the  length  of  the  edge  and  A  is  the  area  of  triangle  T“ .  (Note 
n  n  n 

that  we  use  the  convention,  followed  throughout  the  report,  that  subscripts 

refer  to  edges  while  superscripts  refer  to  faces.)  The  new  basis  function  f^ 

is  to  be  used  to  approximately  represent  the  current,  and  we  list  here  some 

of  its  properties  which  make  it  uniquely  suited  to  this  role: 

i)  The  current  has  no  component  normal  to  any  of  the  edges  except 

the  common  edge  (edge  n)  of  T+  and  T  ;  were  this  not  the  case,  the 

n  n 

continuity  equation  (A)  would  demand  the  presence  of  line  charges 
along  these  edges. 


n  -+ 


2A 


+  n 


p  ,  r  in  T 


—  p  ,  r  in  T' 


2A 


(6) 


,  otherwise. 


ii)  The  surface  divergence  of  the  basis  current,  which  is  proportional 


to  the  surface  charge  density,  is 


-T  ,  r  in  T+ 

A+  n 

n 


V  *f  = 
s  n 


n  ,  r  in  T 

—  -  r 


0  ,  otherwise. 


(7) 


+  +  +  + 

where  the  surface  divergence  in  T  is  (±l/p  )  3(p  f  )/9p  .  The 

n  n  n  n  n 

charge  density  is  thus  constant  in  each  triangle,  the  total  charge 
associated  with  the  triangle  pair  T+  and  is  zero,  and  the  basis 
functions  for  the  charge  evidently  have  the  form  of  a  pulse  doub¬ 
let  [15]. 

iii)  The  component  of  current  crossing  the  nth  edge  is  continuous,  and 
hence  no  line  charge  exists  there;  this  may  be  seen  by  Fig.  3 


which  shows  that  the  normal  component  of  along  edge  n  is  just 

+ 

the  height  of  triangle  T  with  edge  n  as  the  base  and  with  the 

n 


height  expressed  as  (2A  )/i 

n  n 


These  factors  are  used  to  normalize 
f  such  that  its  flux  density  normal  to  edge  n  is  unity  in  (6) , 


hence  ensuring  continuity  of  current  normal  to  the  edge. 


iv)  The  moment  of  f  is  given  by  (A  +  A  )favg 
-  n  n  n  n 


where 


(A+  +  A-)faVg  = 
n  n  n 


f  dS 
n 

t+  +  t" 

n  n 


n  ,-c+  -c->. 

T  (Pr,  +  Pn  > 

z  n  n 


=  l  (?C+-  ?C') 
n  n  n 


(8) 


“  C  — 

and  p  is  defined  between  the  free  vertex  and  the  centroid  of 
n 

+  _p+  _c_ 

T  with  p  directed  away  from  the  vertex  and  p  directed  toward 
n  n  n 
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«ci 

the  vertex,  as  shown  in  Fig.  4,  and  r  is  the  vector  from 

n 

+ 

0  to  the  centroid  of  T';  Eq.  (8)  may  be  most  easily  verified  by 

expressing  the  integral  therein  in  terms  of  area  coordinates, 

which  are  discussed  later  in  this  report. 

Except  for  boundary  edges,  a  basis  function  f  is  associated  with  each  edge 

of  the  triangulated  structure.  The  current  on  S  may  be  approximated  in  terms 

of  the  f  as 
n 

N 

j  =  2  i  i„<;>  <9> 

n=l 

where  N  is  the  number  of  edges  not  on  a  surface  boundary.  Since  a  basis 
function  is  associated  with  each  non-boundary  edge,  up  to  three  non-zero 
basis  functions  may  exist  within  each  triangular  face.  At  each  edge,  how¬ 
ever,  only  the  basis  function  associated  with  that  edge  may  have  a  component 
of  current  normal  to  the  edge;  by  (i) ,  all  other  basis  currents  in  that  face 
are  parallel  to  the  edge.  Furthermore,  since  the  normal  component  of  f^  at 
the  n*"*1  edge  is  unity,  each  coefficient  1^  in  (9)  may  be  interpreted  as  the 
normal  component  of  current  density  flowing  past  the  n^  edge.  Because  the 
normal  component  of  current  at  a  surface  boundary  must  vanish  anyway,  we  need 
not  bother  to  define  basis  functions  associated  with  boundary  edges,  and 
hence  (9)  includes  only  contributions  from  non-boundary  edges. 

The  radial  nature  of  the  current  flow  associated  with  each  basis  function 
is  at  first  disconcerting — certainly  for  a  small  triangle  modeling  a  smooth 
section  of  the  scatterer  surface,  one  would  not  expect  the  direction  of  the 
actual  current  to  vary  substantially  within  the  triangle.  Turning  the  ques¬ 
tion  around,  one  might  ask,  "Can  a  superposition  of  the  basis  functions  within 
a  triangle  represent,  say,  a  constant  vector  in  the  triangle?"  That  the 
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oraetry  of  vectors  to  centroids 
sociated  with  an  edge. 


answer  is  affirmative  can  be  seen  with  the  aid  of  Fig.  5,  which  shows  a 
triangle  T  with  the  edges  arbitrarily  labeled  1,  2,  and  3.  With  the  vectors 


p^,  p 2>  and  p^  as  shown,  the  basis  functions  in  T  are  f^  *  (SL^/2A)  p^,  1=1, 
2,3,  where  A  is  the  triangle  area  and  where,  for  simplicity,  we  assume  that 
the  current  reference  directions  are  out  of  the  triangle  for  each  edge.  It 
is  apparent  from  the  definition  of  f_^  and  the  figure  that  the  linear  combina¬ 
tions  ^2^1  ”  ^1^2  an<*  ^3^1  -  ^1^3  are  constant  vectors  for  every  point  r  in¬ 
side  the  triangle  and  are  parallel  to  sides  3  and  2,  respectively.  Since 
the  two  composite  forms  are  linearly  independent  (i.e.,  non-parallel),  a  con¬ 
stant  vector  of  arbitrary  magnitude  and  direction  within  the  triangle  may  be 
synthesized  by  an  appropriate  linear  combination  of  the  two  forms,  as  asserted. 


Testing  Procedure 

The  next  step  in  applying  of  the  method  of  moments  is  to  select  the  test¬ 
ing  procedure.  As  testing  functions,  we  choose  the  expansion  functions  f 
developed  in  the  previous  section.  With  the  symmetric  product  definition 


<?»  i> 


g  dS, 


(10) 


we  test  Eq.  (5)  with  ?m,  yielding 

<Ii,  f  >  =  jw  <A,  f  >  +  <V$,  f  >.  (11) 

m  mm 

By  standard  surface  vector  calculus  formulas  [19],  the  last  term  in  (11)  can 
be  rewritten  as 


<V<J>,  f  > 

m 


$>V  •  f  dS, 

s  m 


(12) 


where  use  has  been  made  of  the  fact  that  none  of  the  f  has  a  component  nor¬ 
mal  to  any  part  of  the  boundary  of  S.  Using  (7) ,  we  next  approximate  the 
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1 


integral  in  (12)  as  follows: 


4>V 


f  dS 
m 


l  [- 

®  L+ 


$  dS  -  - 


4>  dS 


m  m 


m  m 


=  A  f*  (r®*)  -  *  (r®“)  1 
m  [_  m  m  j 


(13) 


where  the  two  averages  of  $  over  triangles  T  and  T  have  been  approximated  by 

mm 

the  corresponding  values  of  $  at  the  centroids  of  the  triangles.  We  similarly 
approximate  the  integration  of  the  vector  potential  and  incident  field  terms 
in  (11): 


'-i' 

E 


vA  , 


,  £  > 
*  m 


2A 

m  m 
rsi.-c+x 


E1 


•  p+  dS  +  — 
m  2A~ 


1/-C-. 


V] 

—  1 

- 

■  •  dSJ 

;T“ 

A 

m 

V. 

TM*(Jf)r  (P" ) +  {*  rol’ 


(14) 


si 


where  the  integrals  are  eliminated  by  approximating  E  and  A  with  their  values 
at  the  centroid  of  each  triangle  and  then  carrying  out  integrations  similar 
to  those  used  to  obtain  (8).  With  (12)  -  (14),  (11)  now  become.' 


»„  [5<;f  >  •  t + 5<c>,  •  t  ] +  v  -  t<jr>] 

■i  [ i1(r'+)  •  ^  1  , 

ml  m  2  m  2  I 


(15) 


which  is  the  equation  to  be  enforced  at  each  triangle  edge,  m*=  1,2 . N. 

Another  interpretation  of  the  testing  procedure  arriving  at  (15)  is  also  pos¬ 
sible.  One  may  integrate  the  vector  component  of  (5)  parallel  to  the  path 

*— C"t*  — c*f*  ^cf  — ^  ^ 

from  the  point  r  to  (r  +  p  /2)  and  thence  to  r  ,  approximating  E  and  A 
m  m  m  m 
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along  each  portion  of  the  path  by  their  respective  values  at  the  triangle 


centroids.  The  resulting  equality,  when  multiplied  by  £,m>  is  Eq.  (15).  Undi 
either  interpretation,  the  purpose  of  the  testing  procedure  is  reduce  the 
differentiability  requirement  on  $  by  integrating  it  first.  The  purpose  of 
approximations  (13)  and  (14)  is  to  remove  all  surface  integrals  of  potential 
quantities ;  were  this  not  done,  a  prohibitively  expensive  two-fold  surface 
integration  would  be  required  to  fill  the  moment  matrix  since  computation  of 
the  potentials  themselves  already  involves  one  surface  integration. 


Evaluation  of  Matrix  Elements 

Substitution  of  the  current  expansion  (9)  into  (15)  yields  an  N  x  N 
system  of  linear  equations  which  may  be  written  in  matrix  form  as 


Z  I  =  V  (16) 

where  Z  =  [Z  ]  is  an  N  *  N  matrix  and  I  =  [I  ]  and  V  =  [V  ]  are  column  vec- 
mn  n  m 

tors  of  length  N.  Elements  of  Z  and  V  are  given  by 


Z  =  i 
mn  m 


+  A 

mn 


(17) 


V  =  l 
ra  m 


E+  •  Pm  +  E~ 
m  —y—  m 


(18) 


where 


-jkR" 


Jr  =  -r-  f  f  (r')  “  t —  ds'  , 
mn  4tt  J  n  ± 


f  v  •  f 


-jkR' 


m 


dS’ 


(19) 


(20) 


and 


For  plane  wave  incidence,  we  set 

E1^)  =  (E  0  +  E  $  )  ejt  ‘  r 

Bo  <p  o 

where  the  propagation  vector  ic  is 


(21) 


(22) 


k  =  (k  sin  8  cos  4>  x  +  sin  0  sin  6  y  +  cos  6  z)  (23) 

o  o  oo  o 

and  (0  ,  <p  )  defines  the  angle  from  which  the  plane  wave  arrives  in  the  usual 
o  o 

spherical  coordinate  convention.  The  unit  vectors  §  and  $  are  constant 

o  o  - 

vectors  which  coincide  with  the  spherical  coordinate  unit  vectors  at  points 

on  the  line  from  0  in  the  direction  of  k.  Once  the  matrices  Z  and  V  of  (16) 

are  determined,  one  may  solve  the  system  of  linear  equations  for  I. 

We  note  that  although  a  general  matrix  element  Z  is  associated  with 

mn 

the  pair  of  edges  m  and  n,  each  computed  integral  appearing  in  Z  is  actu- 

mn 

ally  related  to  a  source  triangle  attached  to  edge  n  with  an  observation 
point  at  the  centroid  of  a  triangle  attached  to  edge  m.  For  each  such  ob¬ 
servation  and  source  triangle  pair,  these  same  integrals  contribute  to  an 
element  of  Z  whose  row  index  corresponds  to  one  of  the  observation  triangle 
edges  and  whose  column  index  corresponds  to  one  of  the  source  triangle  edges. 

Thus,  rather  than  individually  compute  each  element  of  Z  ,  we  instead  com- 

‘  mn 

pute  all  vector  and  scalar  potentials  associated  with  each  observation-  and 
source- face  combination  and  then  place  the  quantities,  appropriately 
weighted,  into  the  elements  of  Z  corresponding  to  the  various  edges  associ¬ 
ated  with  these  faces.  (The  face  matrix  described  in  Appendix  A  provides  a 
convenient  means  for  keeping  track  of  the  correspondence  between  faces  and 
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edges  as  well  as  for  determining  the  current  reference  directions  within  each 
patch.)  Doing  the  computations  in  this  fashion  results  in  up  to  a  nine-fold 
increase  in  efficiency  in  filling  the  matrix  Z  over  the  direct  edge-by-edge 
approach. 

In  accordance  with  the  above  discussion,  let  us  consider  the  evaluation 
of  the  vector  and  scalar  potential  integrals  for  a  given  source  and  observa¬ 
tion  face  combination.  Fig.  6  illustrates  an  observation  point  in  face  p 
with  current  sources  residing  in  face  q.  For  purposes  of  illustration,  we 
assume  the  edges  of  face  q  are  numbered  1,  2,  and  3  with  edge  lengths  £^, 
and  and  opposite  vertices  at  r^,  r^,  and  r^,  respectively.  We  further 
denote  face  q  simply  as  triangle  T^with  area  A^.  Each  of  the  three  basis 
functions  which  may  exist  simultaneously  in  T  is  proportional  to  one  of  the 
vectors  p^,  p^,  and  p^  defined  in  Fig.  6,  where  the  subscripts  correspond  to 
the  associated  edges  and  we  have  dropped  the  ±  superscripts.  Each  of  the 
vectors  p.. ,  i  =  1,2,3,  is  shown  directed  away  from  its  associated  vertex  in 
the  figure,  but  may  instead  be  directed  toward  the  vertex  if  the  current 
reference  direction  for  that  edge  is  into  the  triangle.  Consequently, 

P±  =  ±  (r'  -  r±),  i  =  1,2,3,  (24) 

where  the  positive  sign  is  used  if  the  positive  current  reference  direction 
is  out  of  and  the  negative  sign  is  used  otherwise.  We  wish  to  evaluate 
the  magnetic  vector  potential. 


and  electric  scalar  potential 


(26) 


associated  with  the  i  basis  function  on  face  q  observed  at  the  centroid 
of  face  p.  In  (25)  and  (26), 


(27) 


—  cp 

where  r  r  is  the  position  vector  of  the  centroid  of  face  p. 

Integrals  (25)  and  (26)  are  most  conveniently  evaluated  by  transformin'g 
to  area  coordinates  [20]  within  the  source  triangle.  Fig.  7  shows  the  posi¬ 
tion  vector  r1  at  some  arbitrary  point  in  Tq.  The  vectors  then  divide  Tq 
into  three  regions  of  areas  A^,  A^ ,  and  A^  which  are  constrained  to  satisfy 
A^  +  A^  +  A^  =  A^.  We  define  the  normalized  area  coordinates  as 


£ 

n 

S 


(28a) 

(28b) 

(28c) 


which,  because  of  the  area  constraint,  must  satisfy 

C  +  H  +  C  =  1.  (29) 

Note  that  all  three  coordinates  vary  between  zero  and  unity  in  and  that  at 
the  triangle  corners,  r^,  r^,  and  r^, the  triplet  (£,  n,  C)  takes  on  the  values 
(1,0,0),  (0,1,0),  and  (0,0,1),  respectively.  The  transformation  from  Cartesian 
to  area  coordinates  may  be  written  in  vector  form  as 

r’  =  ^  +  nr 2  +  Cr3,  (30) 


21 


iPq  _  tpq  _  xpq  _  Tpq  .  (34d) 

Thus  we  see  that  only  three  independent  integrals,  (34a)  -  (34c),  must  be 

numerically  evaluated  for  each  combination  of  face  pairs  p  and  q.  The  three 

integrals,  in  turn,  contribute  to  up  to  nine  elements  of  Z  in  (17).  For  a 

closed  object,  the  number  of  independent  integrals  to  be  computed  turns  out 
2 

to  be  (4/3)  N  .  Numerical  evaluation  of  the  integrals  (34a)  -  (34c)  may  be 
accomplished  by  using  numerical  quadrature  techniques  specially  developed  for 
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triangular  domains  [21]  together  with  the  procedures  discussed  in  Appendix  B. 


III.  MAGNETIC  FIELD  FORMULATION 

In  this  section,  the  magnetic  field  integral  equation  (MFIE)  is  derived 
for  a  conducting  scatterer  S.  Since  the  MFIE  applies  only  to  closed  bodies, 
throughout  this  section  we  assume  that  S  has  no  boundary  edges.  The  vector 
basis  functions  f  of  the  previous  section,  used  there  as  expansion  and  test¬ 
ing  functions  for  the  EFIE,  are  chosen  to  play  the  same  roles  here  in  the  nu¬ 
merical  solution  of  the  MFIE.  The  resulting  moment  matrix  elements  are  given 
and  their  numerical  evaluation  is  also  discussed  in  this  section. 

Magnetic  Field  Integral  Equation 

The  magnetic  field  integral  equation  is  derived  by  noting  that  the  in¬ 
duced  current  3  on  S  is  related  to  the  incident  and  scattered  magnetic  fields 
“■  i  —  g 

H  and  H  ,  respectively,  by 

3  =  n  x  (H*  +  HS),  (35) 

where  n  is  an  outward  unit  normal  vector  on  S.  It  may  be  shown  by  a  detailed 
limiting  argument  [19]  that  for  observation  points  r  not  on  an  edge, 

n  x  hs  =  lim  n  x  V  x  A  . 
r+S 

=  f  +  n  x  J-  j  jx  V'GdS’,  (36) 

where  G  =  exp  (-jkR)/R  and  r  approaches  S  from  the  exterior.  Combining  (35) 
and  (36),  we  obtain  the  magnetic  field  integral  equation  (MFIE): 
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y+ 


*  -i 

nx  H 


J 

2  -n 


Air 


J  x  V'GdS' 


(37) 


Eq.  (37)  is  an  integral  equation  of  the  second  kind  (i.e.,  the  unknown  J  ap¬ 
pears  outside  as  well  as  under  the  Integral),  and  the  kernel  is  regular.  In 
fact,  if  S  is  a  triangulated  surface  and  r  is  some  point  on  S  interior  to  a 
planar  triangular  face,  then  the  current  in  that  face  does  not  contribute  to 

_  _  _  A 

the  integral  since  J  x  (r  -  r')  is  parallel  to  n  there.  A  slight  modifica¬ 
tion  to  (36)  and  (37)  is  required  if  r  is  directly  on  an  edge  [22];  this  sit¬ 
uation  will  not  arise  in  the  present  approach,  however. 


Expansion  and  Testing  Procedure 

As  with  the  EFIE,  we  find  the  functions  F  to  be  suitable  both  as  expan¬ 
sion  and  testing  functions.  Treating  the  testing  procedure  first,  we  test 
(37)  with  f,  and  use  approximation,  paralleling  thoae  yielding  (13)  and  (14) 
to  obtain 
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(38) 


where  n  is  the  outward  surface  normal  in  triangle  T  and 
m  m 

for  J  results  in  a  matrix  equation  of  the  form 

8  I  *  I1, 


Substituting  expansion  (9) 


(39) 
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where  the  elements  of  the  matrix  6  ■  [{3  ]  and  the  column  vector  I*  ■  [I1] 

mn  m 

are  given  by 


6  =  h  <?,  f  >  -  S.  -f-  *  n*  *  TZ 

ran  in  n  in  z  in  htt 


f  x  (V'G)+  dS’ 
n  m 


^m  .-1 
+  — -z-  •  n  x  7— 
2  m  4TT 


f  x  (V 
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•G)-m  dS] 


(40) 


and 


c 


-c- 

p. 


I1  =  l  f  n+  x  Hi(rc+)  •  —  +  n"  X  Hi(rc~)  •  -f-1 .  (41) 

m  m  I  m  m  L  m  m  i 


Solution  of  (39)  yields  the  column  vector  I  of  coefficients  of  the  current 
expansion. 


Evaluation  of  Matrix  Elements 

One  notes  that  a  matrix  element  8  is  associated  with  edge  pair  m  and 

mn 

n  of  S,  whereas  the  integrals  and  observation  points  appearing  on  8^ 
in  (40)  are  associated  with  the  faces  that  are  connected  to  edges  m  and  n. 

Consequently,  it  is  efficient,  as  with  the  EFIE,  to  evaluate  all  integrals 
required  for  a  given  face-face  combination  and  then  to  sort  the  integrals 
into  the  appropriate  rows  and  columns  of  8  with  the  aid  of  the  face  matrix 
(Appendix  A),  which  provides  a  mapping  from  faces  to  edges. 

Referring  to  Fig.  6  and  the  analysis  following  (24)  ,  we  may  write  the 
required  integrals  in  terms  of  the  vector  integral 


-jkRP 

p.  x  (JCP  -  r ' )  (l+jkRp)  dS’,  P*q,  (42) 

1  (RP) 
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which  represents  the  contribution  at  the  centroid  of  the  p*"*1  triangle  due  to 
the  ith  basis  function  (i  =  1,2,3)  in  triangle  q.  If  p=q,  Ipq  =  0.  Eq.  (42) 
can  be  written  in  terms  of  area  coordinates  defined  on  Tq  as 


!r  *  4^ ((?cp  *  ri}  jpq  +  <?i  - ?cp)  *  (?ijpq  +  vsq +  vpq)]’ (43) 


where 
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and  R  is  given  by  (27) .  The  evaluation  of  (44a)  -  (44c)  may  be  accomplished 
by  numerical  quadrature  by  the  method  of  [21], 

We  next  consider  the  evaluation  of  contributions  from  each  patch  to  the 
symmetric  product  term  in  (40) .  If  edges  m  and  n  do  not  lie  on  a  common  tri¬ 
angle,  there  is  no  contribution  to  the  symmetric  product.  If  they  lie  on  a 
common  triangle  Tq ,  then  let  us  assume  for  illustration  that  m,n  =  1,2,  or  3 
as  in  Fig.  6.  Then  the  contribution  from  Tq  to  the  symmetric  product  is 
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+  r  •  r  1  (45) 

m  n  _[ 

where  the  positive  sign  is  selected  if  the  reference  directions  of  edges  m 
and  n  are  both  directed  either  into  or  out  of  Tq;  the  negative  sign  is  se¬ 
lected  otherwise.  If  m=n,  two  terms  of  the  form  of  (45)  contribute  to  8  , 

DOT 

one  from  T+  and  one  from  T  . 

m  m 

IV.  NUMERICAL  RESULTS 

In  this  section,  numerical  results  are  presented  for  current  distributions 
induced  on  selected  scatterers  under  plane  wave  illumination.  The  geometries 
considered  are  a  conducting  square  plate,  a  bent  plate,  a  circular  disk,  and 
a  sphere.  The  first  three  of  these  involve  open  surfaces  and  therefore  test 
the  EFIE  approach  when  edges  are  present.  The  disk  is  also  an  example  of  a 
structure  whose  curved  boundary  is  not  amenable  to  modeling  by  rectangular 
patches,  and  the  sphere  is  both  an  example  of  a  closed  surface,  to  which  both 
the  EFIE  and  MFIE  apply,  and  of  a  doubly-curved  surface,  which  is  not  amenable 
to  rectangular  patch  modeling 

Flat  Plate 

Fig.  8  and  9  show  the  dominant  component  current  distributions  along  the 
two  principal  cuts  on  a  square  plate  illuminated  by  a  normally  incident  plane 
wave.  For  comparison,  the  solution  of  Glisson  [18],  obtained  using  rectang- 
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Fig.  8.  Distribution  of  dominant  component  of  current 
on  0.15A  square  flat  plate. 


Fig.  9.  Distribution  of  dominant  component  of  current 
on  a  1.0A  square  flat  plate. 
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ular  patches,  is  also  given.  The  number  of  patches  listed  in  the  figures  re¬ 
fers  to  the  number  of  charge  patches  in  the  earlier  solution  of  Glisson  and 
to  the  number  of  triangles  (also  equal  to  the  number  of  charge  patches)  in 
the  present  solution.  Note  that  these  quantities  play  similar  roles  in  the 
two  approaches.  No  comparison  of  the  rate  of  convergence  of  the  two  approaches 
should  be  inferred  from  the  figures  since  both  solutions  are  already 
well-converged  for  the  number  oi  unknowns  used.  Note  also  that  the  density  of 
data  points  appearing  in  the  figures  for  the  triangular  patch  solution  is  not 
truly  indicative  of  the  linear  density  of  the  subdomains.  This  is  because, 
in  effect,  we  show  data  points  for  every  other  edge,  i.e.  only  for  those  edges 
where  the  associated  current  normal  to  the  edge  is  parallel  to  the  current 
component  we  wish  to  observe. 

Fig.  8  shows  the  current  induced  on  a  plate  0.15X  on  each  side.  At  this 
low  frequency,  the  current  distribution  is  largely  determined  by  the  edge 
conditions  and  this  case  therefore  provides  a  good  test  of  the  method' s  ability 
to  handle  surface  edges.  We  note  the  absence  of  any  anomalies  in  the  com¬ 
puted  distribution  near  the  plate  edges.  The  elimination  of  such  anomalies 
is  attributed  to  using  basis  functions  ii  which  the  expansion  coefficients 
are  not  associated  with  currents  at  plate  edges  and  to  a  testing  procedure  in 
which  potentials  are  not  evaluated  at  edges  [15]. 

Fig.  9  shows  corresponding  results  for  a  1.0A  square  plate.  It  also 
shows  that  the  edge  behavior  of  the  current  distribution  is  confined  to  a 
smaller  region  near  the  edges  than  for  the  0.15X  plate  and  that  the  current 
on  the  interior  portion  of  the  plate  is  beginning  to  exhibit  the  physical 
optics-nlus-standing  wave  distribution  characteristic  of  the  higher  frequencies. 
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Bent  Plate 


Figs. 10  and  11  show  the  dominant  component  current  distributions  along 
the  center  and  perpendicular  to  the  direction  of  current  flow  at  two  frequencies 
on  a  square  plate  bent  through  an  angle  of  50° •  The  bend  is  located  at  a 
distance  of  one-third  the  plate  width  from  an  edge  and  a  plane  wave  with  the 
electric  field  polarized  parallel  to  the  bend  is  incident  normal  to  the 
larger  section  of  the  bent  plate.  Other  polarizations  and  angles  of  incidence 
have  been  examined  and  the  resulting  current  distributions  show  good  corre¬ 
spondence  with  those  of  Glisson  [18]. 


Circular  Disk 

Fig.  12  shows  the  computed  current  disbribution  on  a  circular  disk  il¬ 
luminated  by  a  normally  incident  plane  wave.  The  component  is  shown  along 
the  cut  across  the  diameter  oriented  perpendicular  to  the  incident  electric 
field  vector.  Also  shown  for  comparison  is  the  quasi-static  solution  valid 
at  low  frequencies  [23]. 


Sphere 

Fig.  13  shows  the  current  distribution  computed  by  the  EFIE  along  the 
principal  cuts  on  a  0.2A  radius  conducting  sphere.  The  cases  of  axial  inci¬ 
dence  and  equatorial  incidence  are  both  considered  in  order  to  observe  the 
influence  of  the  triangulation  scheme  on  the  solution.  Also  shown  for  com¬ 
parison  is  the  exact  eigenfunction  solution.  Both  illuminations  show  very 
good  agreement  with  the  exact  solution. 

Since  the  sphere  is  a  closed  body,  this  problem  may  also  be  examined 
using  the  MFIE.  Fig.  14  shows  the  results  of  the  MF1E  computation.  Compar¬ 
ison  of  the  exact  solution  with  the  MFIE  computation  is  disappointing,  partic- 
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_ 9  6  RECTANGULAR  PATCHES  (GLISSON) 

OO  72  TRIANGULAR  PATCHES 


Fig.  10.  Distribution  of  dominant  component  of 
current  on  a  0.15A  bent  square  plate. 
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Distribution  of  current  on  a  0.05A  radius 
circular  disk. 


6  (DEGREES) 


Fig.  13,  Distribution  of  current  components  on  a 
0.2A  radius  conducting  sphere  calculated 
by  an  electric  field  integral  equation 
formulation. 


36 


Fig.  14.  Distribution  of  current  components  on  a  0.2A 
radius  conducting  sphere  calculated  by  a  mag¬ 
netic  field  integral  equation  formulation. 
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ularly  when  compared  with  the  EFIE  solution.  In  an  attempt  to  try  to  improve 
the  accuracy  of  the  MFIE,  an  alternative  testing  procedure  was  also  examined 
in  which  point-matching  at  the  mid-points  of  edges  was  used.  Although  this 
had  the  effect  of  radically  changing  several  of  the  elements;  including  the 
diagonal.,of  the  matrix,  the  resultant  current  distributions  were  virtually 
indistinguishable  from  Fig.  14  .  Due  to  computer  budget  limitations,  further 
experimentation  with  frequency,  number  of  unknowns,  and  triangulation  schemes 
was  not  possible.  We  point  out,  however,  that  the  surface  discretization 
used  results  in  a  rather  crude  approximation  to  the  sphere  and  we  suggest  that 
perhaps  the  good  agreement  in  the  EFIE  case  may  have  been  largely  fortuitous. 
More  experimentation  is  obviously  needed  to  establish  the  superiority  of 
either  formulation  for  closed  bodies. 


V.  SUMMARY 

In  this  report, the  electric  field  integral  equation  (EFIE)  is  used  with 
the  moment  method  to  develop  a  simple  and  efficient  numerical  procedure  for 
treating  problems  of  scattering  by  arbitrarily-shaped  objects.  The  objects  are 
modeled  for  numerical  purposes  by  planar  triangular  surface  patch  models. 
Because  the  EFIE  formulation  is  used,  the  procedure  is  applicable  to  both 
open  and  closed  bodies.  Crucial  to  the  formulation  is  the  development  of  a 
set  of  special  subdomain  basis  functions  defined  on  pairs  of  adjacent  triangu¬ 
lar  patches.  The  basis  functions  yield  a  current  representation  which  is  free 
of  line  or  point  charges  at  subdomain  boundaries. 

A  second  approach  using  the  magnetic  field  integral  equation  (MFIE) 
and  employing  the  same  basis  functions  is  also  developed.  Although  the  MFIE 
applies  only  to  closed  bodies,  the  moment  matrix  of  the  MFIE  is  also  needed 
in  dielectric  scattering  problems  and  in  the  so-called  combined  field  integral 
equation  used  to  eliminate  difficulties  with  internal  resonances  present  in 
the  MFIE  and  EFIE  formulations. 

The  EFIE  approach  is  applied  to  the  scattering  problems  of  plane  wave 
illumination  of  a  flat  square  plate,  a  bent  square  plate,  a  circular  disk, 
and  a  sphere.  Comparisons  of  surface  current  densities  are  made  with  previous 
computations  or  exact  formulations  and  good  agreement  is  obtained  in  each 
case.  The  MFIE  approach  is  also  applied  to  the  sphere  and  reasonable  agree¬ 
ment  with  exact  calculations  is  obtained. 
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APPENDIX  A 


TOPOLOGICAL  PROPERTIES  AND  MATHEMATICAL 
REPRESENTATION  OF  TRIANGULATED  SURFACES 


In  this  appendix  we  consider  some  topological  properties  of  a  triangu¬ 
lated  surface  and  present  a  simple  mathematical  representation  for  such  a 
surface.  In  the  topological  discussion  a  number  of  geometrical  quantities 
are  defined  and  some  relationships  between  them  are  given.  Consideration 
is  then  given  to  a  means  for  mathematically  representing  a  triangulated 
surface  in  a  form  that  is  convenient  whether  the  surface  data  is  supplied  by 
the  modeler  or  is  generated  by  an  automatic  surface  triangulation  computer 
subprogram  [ 24] .  From  this  representation  may  be  derived  an  alternative 
representation  which  is  actually  more  convenient  for  the  subsequent 
numerical  processing  necessary  in  applying  the  moment  method. 

We  mention  here  at  the  outset  that  one  aspect  of  the  electrical  repre¬ 
sentation  of  the  scattering  problem  also  has  a  bearing  on  what  information 
is  required  in  the  geometrical  representation  of  the  surface.  That  factor 
is  that  the  unknowns  to  be  solved  for  are  the  components  of  current  normal 
to  each  triangle  edge.  There  are  two  possible  senses  in  which  the  current 
can  flow  normal  to  each  edge  and  the  modeler  should  be  able  to  select  either 
choice  to  establish  an  assumed  reference  direction  for  the  current  at  each 
edge.  Furthermore,  it  is  desirable  that  the  reference  current  specification 
be  incorporated,  if  possible,  in  the  geometrical  representation  so  as  to 
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minimize  the  amount  of  user-supplied  data.  As  will  be  seen,  the  geometrical 
representation  given  here  does  indeed  incorporate  the  specification  of  cur¬ 
rent  reference  directions. 


An  arbitrary  body  modeled  by  triangular  patches  is  shown  in  Fig.  Al. 

The  body  is  assumed  to  be  connected,  orientable,  of  finite  extent,  and 
composed  of  non- intersecting  surfaces.  In  general,  a  triangulated  surface 
modeling  an  arbitrary  body  consists  of  planar  triangular  faces , 
vertices ,  and  edges .  These  geometrical  elements  are  illustrated  in  Fig. 
Al. 

An  arbitrary  surface  may  also  have  handles .  Roughly  speaking,  a 

handle  is  a  portion  of  a  surface  which,  if  detached  from  the  surface,  would 

resemble  a  torus  (Fig.  Al) .  Any  closed,  orientable  surface  with  N  handles 

h 

may  be  continuously  deformed,  by  twisting  and  stretching,  into  a  sphere  with 

N,  handles.  These  spheres  with  N,  handles  may  be  thought  of  as  canonical 
h  h 

objects  used  to  classify  all  closed,  orientable  surfaces.  In  deforming 

surfaces  into  spheres  with  handles,  edges  are  permitted  to  pass  through  one 

another,  but  they  may  not  be  broken  or  disconnected.  A  surface  with  no 

handles  is  simply-connected;  a  surface  with  N.  handles  is  said  to  be  (N. +1)- 
-  h  h _ 

connected. 

If  a  surface  is  open,  it  is  bounded  by  one  or  more  boundary  curves 
(Fig.  Al) ,  each  of  which  is  assumed  to  form  a  closed,  non-intersecting 
curve  on  the  surface.  We  may  associate  with  each  boundary  curve  an  aperture, 
which  is  any  simply-connected  surface  having  one  and  only  one  boundary  curve 
congruent  to  the  associated  boundary  curve  on  the  triangulated  surface. 


APERTURE 


BOUNDARY 

CURVE 


>► 


Intuitively,  an  aperture  surface  is  merely  a  surface  which  can  be  used  to 

cover  a  "hole"  in  another  surface.  (Unfortunately,  the  term  "hole"  is  not 

appropriate  to  describe  regions  such  as  that  exterior  to  the  boundary  curve 

of  a  rectangular  plate,  and  hence  we  use  the  term  "aperture"  to  describe  a 

closing  surface  instead.  For  the  rectangular  plate,  for  example,  a  suitable 

aperture  is  a  rectangular  box  with  one  open  end.)  We  assume  that  an 

arbitrary  surface  has  apertures  (and  associated  boundary  curves)  and 

that  a  total  of  edges,  called  boundary  edges,  lie  on  these  boundary 

curves.  If  there  are  no  boundary  edges,  N  =0  and  the  surface  is  closed. 

b 

We  next  employ  a  theorem  due  to  Euler  which  states  that  [25]  for  closed 
surfaces , 

N'  -  N'  +  N'  =  2(1-0  .  (Al) 

rev  h 

The  primes  remind  us  that  the  result  applies  only  to  closed  surfaces.  The 
right  hand  side  of  (Al)  is  a  topological  invariant  known  as  the  Euler 
characteristic  and  is  the  same  for  any  closed  surface  which  can  be  contin¬ 
uously  deformed  into  a  sphere  with  handles.  To  extend  the  theorem  to 
open  surfaces,  we  first  close  all  the  apertures.  This  may  be  accomplished 
by  introducing  for  each  aperture  an  auxiliary  vertex  and  auxiliary  edges 
connected  between  this  vertex  and  each  vertex  on  the  associated  boundary 
contour  (Fig.  A2).  These  auxiliary  vertices  may  be  arbitrarily  located, 
provided  they  do  not  coincide  with  one  another  or  rest  on  edges  or  vertices 
of  the  original  surface.  The  resulting  closed  triangulated  surface  consists 

of  N'  =  N,  +  N,  faces,  N'  =  N  +  N,  edges,  and  N'  =  N  +  N  vertices, 
ffb  eeb  vva 

Substituting  these  relationships  into  (Al),  one  obtains 
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AUXILIARY 

VERTICES 


Fig.  A2.  Auxiliary  vertices  and  edges  used  to  close 
apertures  to  form  a  closed  body. 


N,  -  N  +  N  =  2 (1-N  )  -  N  , 
f  e  v  ha 


(A2) 


an  extension  of  (Al)  to  open  bodies.  We  may  also  eliminate  either  N-  or  N 

f  e 

in  (A2)  by  noting  that  since  all  of  the  faces  are  triangular,  3N'  =  2N'. 

f  e 

This  follows  from  the  fact  that  for  a  closed  surface,  each  edge  is  counted 
twice  if  we  sum  the  three  edges  per  triangle  over  all  the  triangles.  Use 
of  this  relation  to  eliminate  N^.  and  N^,  respectively,  in  Eq.  (A2)  yields 


N  =  3N  +  3N  +  6(N  -1)  -  N, 
e  v  a  h  b 


(A3) 


and 


N  =  2N  +  2N  +  4(N,  -1)  -  N  . 
f  v  a  h  b 


(A4) 


Eq.  (A3)  may  be  used  to  determine  the  number  of  unknowns  to  be  found  if  one 

knows  the  number  of  vertices,  apertures,  handles,  and  boundary  edges  of  the 

model.  Since  currents  normal  to  boundary  edges  are  zero  and  hence  are  not 

solved  for,  the  number  of  unknowns,  N,  is  equal  to  the  number  of  surface 

interior  edges,  N  -  N,  . 

e  b 

A  computer  subroutine,  GEOM,  has  been  developed  to  accept  data  describ¬ 
ing  a  triangulated  surface  and  to  generate  auxiliary  information  and  data 
necessary  for  further  numerical  processing.  The  subroutine  requires  two 
sets  of  input  data.  The  first  is  an  indexed  list  or  vertex  matrix  of 
position  vectors  r  =  (x^,  y^,  z  ^) ,.  i  *  1,  2,  ••*,  Nv>  The  components  of 
the  vectors  r^are  the  Cartesian  coordinates  of  the  itfl  vertex  with  respect 
to  a  global  coordinate  system.  The  second  set  of  data  is  an  edge  matrix 


E  -  [e  ],  i  =  1,  2,  •••,  Ne;  j  *  1,  2 
,  th 


(A5) 


in  which  is  listed  in  the  i  row  the  indices  of  the  two  vertices  to  which 
the  ith  edge  is  connected.  The  order  of  appearance  of  the  vertex  indices 
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in  each  row  of  E  assigns  a  vector  direction  to  each  edge,  with  the  first 
index  corresponding  to  the  tail  and  the  second,  to  the  head  of  the  vector. 
The  cross-product  of  this  vector  with  the  surface  normal  of  a  face  adjacent 
to  the  edge  then  gives  a  positive  reference  direction  for  the  surface 
current  in  that  face.  GEOM  does  not  actually  compute  surface  normals  to 
determine  the  reference  directions,  but  uses  a  procedure  to  be  discussed 
later . 

The  vertex  and  edge  matrices  together  completely  determine  the  surface 
geometry,  the  interconnections  of  the  edges  to  form  triangles,  and  the 
current  reference  directions.  However,  one  notes  that  in  filling  the  moment 
matrix,  one  integrates  over  the  surface  faces  and  that  the  results  of  the 
integrations  must  be  placed  in  rows  and  columns  of  the  moment  matrix 
corresponding  to  the  appropriate  edges.  Hence  it  is  convenient  to  introduce 
a  face  matrix, 

F  =  [ f  1 ,  i  =  1,  2,  •••,  Nf ;  j  «  1,  2,  3,  (A6) 

relating  edges  to  the  corresponding  faces.  The  ith  row  of  the  face  matrix 
contains  the  edge  numbers  of  those  edges  making  up  the  boundary  of  the  i*''1 
face.  Subroutine  GEOM  makes  use  of  information  in  the  edge  matrix  to 
find  each  face,  to  assign  it  an  index,  and  to  fill  out  the  corresponding 
row  of  the  face  matrix.  The  order  in  which  the  numbered  edges  appear  in 
each  row  may  be  used  to  assign  an  orientation  to  the  boundary  curve  formed 
by  the  edges  of  each  face.  We  may  associate  this  orientation  ,  in  turn, 
with  the  direction  of  the  surface  normal  of  each  face  through  the  usual 
convention  relating  a  surface  normal  to  an  oriented  boundary  contour. 

If  the  surface  normals  for  all  the  faces  are  to  be  on  the  same  side  of  the 
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surface  (which  is  always  possible  if  the  surface  is  orientable) ,  then  if 
one  travels  in  the  prescribed  direction  around  the  boundaries  of  two  adjacent 
faces  he  must  traverse  their  common  edge  in  opposite  directions  (Fig.  A3). 
GEOM  makes  use  of  this  property  to  correctly  order  the  elements  in  the  face 
matrix  so  that  the  orientation  of  all  the  face  normals  is  toward  the  same 
side  of  the  surface.  The  direction  of  the  surface  normal  is  itially 
chosen  by  locating  the  lowest  numbered  edge  connected  to  edge  i  =  1. 

Then  these  two  edges  are  treated  as  vectors  pointing  away  from  their  common 
vertex  and  their  cross-product  is  computed  with  edge  i  =  1  as  the  second 
vector  in  the  product.  The  surface  normal  is  then  assumed  to  be  parallel  to 
this  cross-product.  Thus,  by  properly  numbering  the  edges  connected  to 
edge  i  =  1,  the  modeler  can  fix  the  choice  of  the  normal  for  an  open  surface. 
Furthermore,  as  already  noted,  by  properly  ordering  the  elements  of  the 
edge  matrix,  he  may  also  choose  the  positive  reference  direction  for  each 
individual  edge. 

If  the  surface  is  closed  (N  =  0) ,  the  normal  pointing  into  the  region 

b 

exterior  to  the  surface  should  always  be  chosen.  In  this  case,  GEOM 
automatically  determines  whether  the  correct  choice  has  been  made  by  cal¬ 


culating  the  volume  of  the  model  according  to 

V  =  ]  dV  =  (  V* (xx)  dV 


i=i 
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N , 


=  7  n1  x1  a\ 

i-i  x  C 


(A7) 


where  T1  is  the  triangular  patch  formed  by  the  i*"^1  face,  is  the  x- 

component  of  the  normal  to  T'*',  x'*’  is  the  x-coordinate  of  the  centroid  of  T 

c 

with  respect  to  the  global  coordinate  system,  and  A1  is  the  area  of  T  .  If 
the  volume  computed  from  Eq.  (A7)  turns  out  to  be  negative,  then  the  modeler 
has  erroneously  chosen  the  interior  normal  and  GEOM  rectifies  this  by  inter¬ 
changing  the  first  two  elements  in  every  row  of  the  face  matrix.  The 
orientation  information  is  used  in  subsequent  calculations  in  the  following 
way:  If  the  boundary  of  a  face  and  one  of  its  edges  have  the  same  orien¬ 

tation  (as  determined  by  the  face  and  edge  matrix,  respectively)  ,  then  the 
positive  reference  direction  for  the  current  normal  to  the  edge  is  out  of 
the  face.  Otherwise,  the  reference  direction  is  into  the  face  (see  Fig.  A3). 

In  the  course  of  computing  the  face  matrix,  GEOM  also  determines  all 
boundary  edges.  It  a 1 so  sorts  them  into  their  corresponding  boundary  curves 
and  hence  determines  how  many  apertures  are  present.  Then  using  Eq .  ( \2) , 
it  determines  iiow  many  handles  the  surface  has  and  computes  the  number  of 
interior  edges  ( i . e. ,  the  number  of  unknown  current  coefficients)  by  N  = 

N  -  N  .  Finallv,  as  a  partial  check  on  the  correctness  of  the  triangulation 
e  b 

scheme,  GEOM  checks  to  see  if  (A3)  and  (A4)  are  satisfied.  This  test  would 
fail,  for  example,  if  through  an  error  In  an  entry  in  the  connection  matrix 
a  patch  was  not  triangular,  but  quadrilateral. 
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APPENDIX  B 


EVALUATION  OF  INTEGRALS  APPEARING 
IN  THE  ELECTRIC  FIELD  FORMULATION 


In  the  solution  of  the  EFIE  using  triangular  patches,  the  numerical 
evaluation  of  the  following  three  integrals  is  required: 
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where  R  =  |r  -  r' |  and  r'  is  a  position  vector  within  a  triangle  whose  area 
coordinates  are  £  and  n-  If  the  observation  point  r  is  within  the  triangle, 
then  R  will  be  zero  for  some  value  of  £  and  0  and  the  integrands  of  all  three 
integrals  will  be  singular.  To  circumvent  difficulties  with  this  case,  we 
rewrite  (Bl)  -  (B3)  as 
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and 


Sfr- 


I 

n 


l 

n=o 


- jkR  . 

(e-  p  -1)  d^dn  + 

=o  R 


C 


l 

n=o 


r1-r>  „ 

I  f  d^dn. 
JC=o  R 


(B6) 


The  integrand  of  each  of  the  first  terms  in  (B5)  -  (B6)  is  non-singular 
and  hence  can  be  numerically  integrated  by  using  quadrature  formulas  for 
triangular  regions  obtained  by  Hammer  et  al.  [21].  The  remaining  integrals 
in  (B4)  -  (B6)  are  evaluated  analytically  by  the  following  procedure. 

We  begin  by  expressing  r'  in  terms  of  area  co-ordinates  as 

?'=?!+  +  (r3-?1)n,  (B7) 

where  r^,  r^  and  r^  are  the  position  vectors  of  the  vertices  of  the  source 
triangle.  Hence, 


where 


and 


_  _  1  __  __  __  __ 

R  =  | r-r  |  =  | (r-r1>-(r2-r1)  £  -(r3~r1)  n | 

=  [a£2  +  Bn2  +  c£  +  Dn  +  ECn  +  f]\  (B8) 

A  =  ^2'^l!2, 

B  =  i  r3-r;1|  2  , 

c  -■  -  2(r-rj)  •  (?2-?1),  (B9) 

D  =  -  2(r-^)  •  (?3-?1), 

E  =  2(?2-?1)  •  (?3-?i>. 

F  =  ImJ2. 


With  these  definitions,  the  second  integral  in  (B4)  can  be  expressed  as 
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and 


E2  =  - 


(r2~^2>  * 


lv?il 


The  remaining  integrals  in  (Bll)  are  both  of  the  form 


fl 

1  /  2 

[  £n 

Jn=o 

^  /  AjH  +  B^  +  C1  +  D1n  +  E1 

With  the  subsitution  x  =  D  n  +  E  ,  equation  (B12)  reduces  to 


I 


A 


Dl+El 


x=E, 


in 


+  bx  +  c  +  x 


dx, 


where 


a  =  Al/Dl  * 


b  = 


BiDr2AiEi 


c  = 


C1D12-B1E1D1  +  A1E12 

7} 


and,  in  terms  of  the  variable 


x,  R  =  /  ax2  +  bx  +  c. 


Eq.  (B13)  may 


as 


rVEi 

D.  I  =  I  = 

1  A  1  „ 

x=E, 


£n  (R+x)  dx 


and  hence  the  problem  is  reduced  to  evaluating  the  integral 


'l  " 


itn  (R+x)  dx. 


By  some  elementary  substitutions  and  integration  by  parts,  I  may  be 


(B12) 

(B13) 


be  written 

(B14) 

(B15) 

written 
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where 


U2  “  X2  +  2(a-l)  * 

2 

We  note  at  this  point  that  d  is  always  greater  than  or  equal  to  zero, 

2  2 

as  can  be  easily  verified  by  noting  that  R  -  x  is  greater  than  or  equal  to 
zero  for  all  x.  The  last  two  integrals  in  (B17)  may  be  integrated  with  the 
aid  of  tables  and  the  result  combined  with  (B16)  to  yield 


I  =  [xS,n  (R+x)-R-x-d  tan  ^ 

1  d 


2  2 

_ b _  .  ,R  , 

4 (a-1)  ln  (  a-1  )  ] 


+  I3* 


where 


fU2  uR 
"  Ju,  u2  + 


du. 


With  the  substitution  z  =  u  +  jd,  can  be  written  as 


(B18) 


I 


3 


Re 


Using  Eq.  2.267.1  of  [26]  and  taking  the  real  part  with  some  straight-forward 
but  tedious  algebra,  one  obtains 


{  X  +  2(fcl) _  *n  (R+x) 


-  jin  \2/a  R+2ax+b| 

l/a.  (a-1) 


+  d 


tan 


ix 

d 


tan 


2dR(a-l)  I] 
bx  +  2c  _  / 


(B19) 
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FT 


Combining  the  results  of  all  these  steps,  one  may  now  write  l'as 


where 


'  =  - I -  l  (p1-)-  |u  £n  (R  +  x  ) 

IJ.-5.I  i-1  °i  li  1  1 


2  1 


b. 


— -  £n  |2/a7  R.  +  2a  x  +  b .  | 

2/17  (arl)  11  11  1 


+  tan 


-1 


(B20) 


h*o 


R.  =  /  A^2  +  +  Ct, 


U  =  x,  + 


i  i  2(ai-l)  » 


X.  =  +  Ei  , 


Ai  =  I r 3~r t (  , 


=  -2(r;J-ri)  •  (r-j^). 


Ci  =  f  ?_?i 1 2  ’ 


(VV  *  (?3“V 


2-?ll 


<Vri>  *  (?_?i) 


I  r2-rl' 
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ai  =  Ai/Di2  ’ 

bi  =  <BiDi  ~  2AiEi)/Di2’ 


Ci  =  (Ci°i 


-  BiDiEi+AiEi2>/Di2> 


V1 


4(at-l)2 


Some  comments  concerning  the  evaluation  of  (B20)  are  in  order:  1)  If  either 
or  is  zero  (which  happens  whenever  the  corner  of  the  triangle  at  r  ,  or 
r2>  respectively,  forms  a  right  angle),  Eq.  (B20)  cannot  be  evaluated  in  the 
form  indicated.  A  simple  way  to  circumvent  this  difficulty  is  to  cyclically 
permute  the  assignment  of  vectors  r^,  r2  and  r^  to  the  corners  of  the  triangle 
until  neither  nor  is  zero,  i.e.  until  the  right  angle  corner  is  placed 
at  r^.  Under  this  new  assignment  of  vertex  indices,  (B20)  is  valid.  2)  The 
argument  of  either  of  the  logarithmic  terms  appearing  in  (B20)  may  vanish 
for  certain  combinations  of  parameters.  Whenever  this  situation  occurs,  how¬ 
ever,  the  corresponding  coefficient  of  the  logarithmic  term  also  vanishes  and, 
by  L' Hospital's  rule,  the  product  can  be  shown  to  vanish  as  well. 

At  this  stage,  all  that  remains  is  to  evaluate  the  singular  integrals  in 
(B5)  and  (B6) .  Consider  the  pair  of  integrals. 


and 


XP  = 


l-n 


CdCdU 

R 


n=o  f=o 


(B21a) 
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l-n 


u-o  c=o 


uagdu 

R 


where 


=  j^A?2  +  Bn2  +  C£  +  Dr)  +  E£n  +  fJ^. 


(B21b) 


With  =  A,b^  =  C  +  Eri  and  =  Br)  +  Dri  +  F,  we  have 


r  2 

R-  Lai?  +bi5  +  cij  • 


Using  Eq.  380.011  of  [27],  one  obtains  for  (B21a) , 
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n=o 


n  (A+B-E)  +  n  (-2A-C+D+E)  +  A+C+F 


»s 


du 
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[Bn2  +  Du  +  F]5*  du  -  -~ 
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fl 


l-n 
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d£du 
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u=0  £=0 


Udgdu 
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(B22) 


One  notes  that  the  last  two  integrals  in  (B22)  are  I',  evaluated  earlier,  and 
Ip  of  (B21b).  Again  using  Eq.  380.201  of  [  27]  and  evaluating  the  first  two 
integrals  in  (B22),  one  obtains 


J1“J2 


—  i ' - i 

2A  2A  Q 


(B23) 


where 
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Equations  (B23)  and  (B24)  may  be  solved  simultaneously  for  the  unknown  inte¬ 
grals  Ip  and  I  to  obtain 

4B(J  -J  )  -  2E(J3-J4)  -  (2BC-ED)  I' 

Ip  2 

4AB-E 

and 

4A(J  -J.)  -  2E(J  -J  )  -  (2AD-EC)  I' 

T  _  3  4 _ 1  2 _ 

Q  (4AB-E2) 

This  completes  the  analytical  evaluation  of  the  singular  integrals. 


(B25) 


(B26) 
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